What is MGcount

MGcount is a RNA-seq quantification tool conceived to incorporate ambiguous read alignments in a flexible way that is compatible with any biotype. This allows to extract more information from total RNA-seq datasets by simultaneous quantifying coding and non-coding transcripts, both small and long, where dealing with heterogeneous read-to-feature alignment ambiguities is key. When aligning reads to a reference genome, we distinguish between two types of ambiguous alignments:

• Multi-mappers: reads aligning to multiple genomic locations.
• Multi-overlaps: reads aligning to a genomic location with multiple annotated features.

To deal with the most frequent multi-overlaps (occuring when a smallRNA loci is embeded within a longRNA), MGcount hierarchically assigns reads to smallRNA and longRNA accounting for their transcript length disparity. Subsequently, MGcount models read-to-feature alignments in a graph. This is exploited to detect and report expression of sequence-similar annotated loci, where reads systhematically multi-map, as integrated features called communities.

System requirements

MGcount is built on the top of FeatureCounts, a well-known computational efficient software (Liao et. al., 2014). Please download it from the following link: http://subread.sourceforge.net/

Installation

MGcount is written in Python and is executed from the command line. You can either download the executable version (single binary file) or install it as a Python3 module.

pip3 install git+https://github.com/hitaandrea/MGcount.git

python3 -m mgcount [args]

chmod +x mgcount

MGcount [args]

1. Hierarchical assignation

First, a hierarchy is used to solve multi-overlapping ambiguities across biotypes during alignment-to-feature assignment based on transcript body length. Reads are assigned to genomic annotated features in three pre-defined sequential rounds: small, long_exon and long_intron. Hence, all reads with at least one alignment overlapping with an annotation in the current round are assigned to such feature and skipped in subsequent rounds.

2. Multi-loci communities recognition

To quantify multi-mapping reads, MGcount builds a directed weighted graph {G=(V,E)} where each vertex (V) is an annotated feature and a pair of directional edges (E) connect two features for which common multi-mapping reads exist. Edge weithgs are defined as the ratio of multi-mapping reads between the two vertex normalized by the total number of reads assigned to the source vertex. Vertex weights are defined as the log-transformed number of assigned alignments (Fig c). Resultant graphs structures capture the multiple-loci topologies of different RNA biotypes. Two graphs are sepparately build for smallRNA and longRNA, using the full pool of input alignments. Subsequently, highly-related features are grouped together by minimizing the map equation with the communities detection approach described by Rosvall et al., 2008 .

3. Expression matrix generation

MGcount generates an output expression matrix for each hierarchical assignation round that are appended together in a single output matrix. For each read, each alignment first gets an 1/N count, where N is the number of multi-mappers or residual multi-overlaps within biotypes that survive the hierarchical assignment. Next, counts for annotations which have been aggregated together in a community by the map equation are summed up. In this way, the systematic ambiguity in multi-mapping reads gets collapsed into a single MG community while the remaining signal is reported as fractionated counts.

Inputs

Three inputs are required to run the program:

• Input alignment file: a .txt file listing the paths to the .bam alignment input files by line
• Annotations gtf file: a .gtf file containing a set of RNA feature annotations
• Output directory path: a string specyfing the path to the directory where MGcount outputs will be stored

Configurable arguments

Outputs

At the end of its execution, MGcount provides the following outputs:

• Count matrix: Expression estimation matrix where each row corresponds to a transcript feature and each column corresponds to 1 input .bam file.
• Features metadata: It in, which match row names in the count matrix; the counting round when assigned; a flag designing whether a feature is a unique annotation or a MG community; and the feature biotype.
• Multi-mapping graph adjacency matrix: Adjacency matrix of links between annotated features allowing to generate the smallRNA or/and longRNA multi-mapping graph, encoded as sparse Matrix in MatrixMarket standard format.
• Multi-mapping graph communities (MGcommunities): Table linking each original feature in the gtf with the resultant feature matching count matrix and feature metadata feature identifiers. It includes both unique features (that remain unmodified) and aggregated features (that are collapsed following MG communities). In addition, the table provides with the total number of alignments per feature.

T1 - Prepare inputs and execute

In the following tutorials, we use two Human Brain samples as example to walk through MGcount execution. First of all, let’s create a folder and download the aligment bam files of the two samples. (The downloading process might take a few minutes).

mkdir mgcount_tutorial
cd mgcount_tutorial

wget https://filedn.com/lTnUWxFTA93JTyX3Hvbdn2h/mgcount/humanbrain_bamfiles.zip
unzip humanbrain_bamfiles.zip -d input_bamfiles

To run MGcount, we need to provide the software with a .txt file specifying the paths to the input aligment files. Here, these are the two .bam alignment files we just downloaded. We can generate this file from the command line:

printf "$PWD/%s\n" input_bamfiles/* > input_bamfilenames.txt

The other required input is a .gtf file with transcript features annotations. MGcount repository provides with four ready gtf files integrating annotations from several databases. Next, we will download them and use the Human annotations file.

wget https://filedn.com/lTnUWxFTA93JTyX3Hvbdn2h/mgcount/integrated_annotations_gtf.zip
unzip integrated_annotations_gtf.zip -d annotations_gtf

Once we have both the gtf and the input .bam files, we can simply run MGcount as a python3 module. A python module is run by calling python3 with the parameter -m and the name of the module. After this, we can specify MGcount arguments.

For this example, we parse the ready-integrateg gtf for the Human genome, the txt file we just created containing the path to the input aligment .bam files and a string designating the directory where MGcount outputs will be generated.

python3 -m mgcount --gtf annotations_gtf/Homo_sapiens.GRCh38.gtf --outdir outputs --bam_infiles input_bamfilenames.txt

If you have been able to run MGcount successfully, your output directory should contain 6 new files generate by MGcount including the RNA count matrix, the feature metadata table and two files containing the graph structure and the communities detected for longRNA and smallRNA respectively.

We are done!

T2 - Explore quantification outputs

In this tutorial, we will load the MGcount outputs we obtained in tutorial 1 and we will explore them from R. MGcount repository contains a few supporting scripts in R providing with side functionalities (manage annotations, visualize MGcount outputs…). It is possible to download the sole folder containing these scripts with the following shell command:

cd mgcount_tutorial
svn export https://github.com/hitaandrea/MGcount/trunk/R

Once the scripts are downloaded, we are ready to launch R. Please, start an R session and define the tutorial root directory as a variable.

root_dir <- '~/mgcount_tutorial/'

Further, this tutorial uses the following R packages. Please, install them with install.packages() if you wish to run this tutorial on your system.

library(dplyr)
library(Hmisc)
library(ggplot2)
library(ggpubr)
library(summarytools)

## Warning in fun(libname, pkgname): couldn't connect to display ":0"

source(paste0(root_dir,'R/integrate_gtf_annotations.R'))

The main output of MGcount is the count_matrix.csv file containing the feature by sample expression matrix. We can import it into R as any regular csv.

counts <- read.csv(paste0(root_dir,'outputs/count_matrix.csv'), row.names = 1)
colnames(counts) <- sub('_Aligned.genome.dedup','',colnames(counts))

Each row in the count table is a feature for which expression has been quantified and each column is associated to a sample. By interrogating for the matrix dimension, we see the matrix contains two columns corresponding to the two human brain libraries. The matrix is ready to be used in any RNA-seq downstream analysis.

dim(counts)

## [1] 47539     2

We can look at the total number of reads assigned from each library by suming up each of the two rows. The mean counts over the two libraries is approximately 4.2 million.

colSums(counts)

## Human_Brain_total_100ng_1 Human_Brain_total_100ng_2 
##                   3837246                   4578696

print(paste('Mean counts:',mean(colSums(counts))))

## [1] "Mean counts: 4207971.05"

Let’s import now the feature_metadata output. This table reports feature-related attributes such as the assignation round to which the annotation belongs, a flag stating wheter the feature is an aggregated community of annotations or an individual feature and the biotype associated to the feature. The feature identifiers in the counts_matrix and the feature_metadata match and therefore, this table can facilitate the extraction of particular features from the count matrix, e,g, a certain biotype, the exonic counts, the subset of features aggregated in communities, etc…)

feat_metadata <- read.csv(paste0(root_dir,'outputs/feature_metadata.csv'))

Here, for example, we profit from the feature_metadata table to look at counts distribution by assignation round.

df <- cbind('counts' = rowSums(counts), feat_metadata)
ggplot(df, aes(x = assignation_round, y = counts)) + geom_violin(fill = 'grey') + 
   scale_y_continuous(trans = 'log', lim = c(10,100000), breaks = c(10,100,1000,10000,100000)) + theme_pubclean() +
  geom_point(position = position_jitter(seed = 1, width = 0.4), size = 0.005, alpha = 0.2)

Also, we can exploit the feature_metadata table to extract stats from the MGcount run. Here we look at the proportion of aggregated features (community_flag variable) by biotype (feature_biotype). The output shows how smallRNA biotypes tend to be more aggregated in communities because duplicated loci are more frequent.

kable(with(feat_metadata, ctable(feature_biotype, community_flag)))

Below, we display a few random rows from the table for illustration.

feat_subset <- with(feat_metadata, feat_metadata[c(
   sample(which(assignation_round == 'small' & community_flag == "True"), 2),
   sample(which(assignation_round == 'small' & community_flag == "False"), 2),
   sample(which(assignation_round == 'long_exon' & community_flag == "True"), 1),
   sample(which(assignation_round == 'long_exon' & community_flag == "False"), 1),
   sample(which(assignation_round == 'long_intron' & community_flag == "True"), 1),
   sample(which(assignation_round == 'long_intron' & community_flag == "False"), 1)),])
kable(feat_subset) %>% scroll_box(height = "300px")

Next we show how to generate a barplot showing the read distribution by biotype group. First lets load a table to group biotypes by category. We will use this to group less abundant biotypes in larger groups for visualization purposes.

bcats <- define_bcats()

feat_df <- merge(feat_metadata, bcats, by.x = 'feature_biotype', by.y = 'biotype', all.x = TRUE, all.y = FALSE) 

## Add exon/intron distinction to biogroup based on counting round
feat_df$biogroup <- as.character(feat_df$biogroup)
feat_df$biogroup[feat_df$biogroup == 'Protein_coding'] <- 'Protein_coding_exon'
feat_df$biogroup[feat_df$biogroup == 'Long_non_coding'] <- 'Long_non_coding_exon'  
feat_df$biogroup[feat_df$assignation_round == 'long_intron'] <-
  sub('exon','intron',feat_df$biogroup[feat_df$assignation_round == 'long_intron'])

feat_df$biogroup <- as.factor(feat_df$biogroup)
feat_df <- feat_df[match(rownames(counts), feat_df$feature),]

We then combine feature_metadata table with the count_matrix again and sum up the counts to get the total number of reads by biotype. We do this sepparately by biotype groups and small-non-coding individual biotypes to further represent the small non-coding spectrum.

## Generate biotype matrix
df <- data.frame(counts %>% group_by(feat_df$biogroup) %>% summarise_all(sum), 
                 check.names = FALSE); names(df)[1] <- 'biotype'
biotype <- reshape(df, idvar = 'biotype', varying = list(names(df)[-1]),
                      v.names = 'counts', times = names(df)[-1],
                      timevar = 'sn', direction = 'long')  
biotype$biotype <- as.character(biotype$biotype)
bgroups <- c("Hybrid","Long_pseudogenes","Protein_coding_intron","Long_non_coding_intron",
             "Protein_coding_exon","Long_non_coding_exon" ,"Short_non_coding","tRNA","rRNA")
biotype$biotype[biotype$biotype %nin% bgroups] <- 'Hybrid'
biotype$biotype <- factor(biotype$biotype, levels = bgroups)

## Generate non-coding biotype table
df <- data.frame(counts[feat_df$biocat == 'sNC',] %>%
                 group_by(feat_df$feature_biotype[feat_df$biocat == 'sNC']) %>% summarise_all(sum), 
                 check.names = FALSE); names(df)[1] <- 'biotype'
biotype_snc <- reshape(df, idvar = 'biotype', varying = list(names(df)[-1]),
               v.names = 'counts', times = names(df)[-1],
               timevar = 'sn', direction = 'long')  

Once we have extracted the counts by biotype group, let’s employ ggplot to visualize the read distribution profiles as barplots.

Reads distribution by biotype:

## ---- Abundance plot by biotype group
colP <- c('violetred1','slateblue1','darkgrey','lightgrey',
          'springgreen4','springgreen3','violetred4','tan3','tan4')
names(colP) <- bgroups
p1 <- ggplot(biotype, aes(x = sn, y = counts, group = biotype, fill = biotype)) +    
  geom_bar(stat = 'identity', colour = 'black', width = 0.8) + 
  coord_flip() + xlab('') + ylab('Number of assigned reads') + theme_pubclean() + 
  guides(fill=guide_legend(nrow=4)) + scale_fill_manual(values=colP) +  
  theme(legend.position = 'top', legend.title = element_blank())

## ---- Relative abundance plot by small non-coding biotype
colP <- c('#e6be97','#848CFF','#4c2382', '#50eb76','tomato','#FDFF87','#FFAE51','darkorchid1','gold1','sienna3')
p2 <- ggplot(biotype_snc, aes(x = sn, y = counts, group = biotype, fill = biotype)) +
  geom_bar(stat = 'identity', position = 'fill', colour = 'black', width = 0.8) + theme_pubclean() +
  coord_flip() + xlab('') + ylab('Proportion of small non-coding assigned reads') + 
  guides(fill=guide_legend(nrow=4)) +   scale_fill_manual(values = colP) +
  theme(axis.text.y = element_blank(), legend.position = 'top', legend.title = element_blank())

## Display plots together
ggarrange(p1, p2, ncol = 2, widths = c(2,1))

Finally, we import the multigraph_communities tables. These tables link each original feature in the gtf with the resultant feature matching the count matrix and feature metadata identifiers. It includes both unique features (that remain unmodified) and aggregated features (that are collapsed following MG communities). Thus, we can track back the features groupped in each aggregated feature.

Here we will use the feature_metadata subset from before to explore the original features in the gtf forming each new MG community feature.

csmall <- read.csv(paste0(root_dir,'outputs/multigraph_communities_small.csv'))
clong <- read.csv(paste0(root_dir,'outputs/multigraph_communities_long_exon.csv'))

kable(subset(csmall, feature %in% feat_subset$feature)) %>% scroll_box(height = "300px")

kable(subset(clong, feature %in% sub('-exon','',feat_subset$feature))) %>% scroll_box(height = "300px")

kable(subset(clong, feature %in% sub('-intron','',feat_subset$feature))) %>% scroll_box(height = "300px")

We are done!

T3- Explore output multi-graph

Below, we use a few functions defined in “mg_visualize.R” to graphically explore the multi-mapping graph topology.

The function ‘mg_build’ takes the Multi-Graph adjacency matrix and the MG communities and creates an igraph object. The next two functions provide with the code to generate a few default plots given a Multi-Graph igraph object.

library(Matrix)
library(igraph)
library(plotly)
source(paste0(root_dir,'R/mg_visualize.R'))

As an example, we will load the smallRNA Multi-Graph adjacency matrix and we will subset it to explore the sub-graph associated to microRNA features.

The adjacency matrix is stored as a sparse symetric matrix in MatrixMarket format, which can be imported in R by the readMM function from ‘Matrix’ R package. We also import the table containing all annotated features and MG communities outputs.

dir.create(paste0(root_dir,'mgplots'))

## Import ml data and adjacency matrix for small assignation round
inputM <- readMM(paste0(root_dir,'outputs/multigraph_matrix_small.mtx'))
multiloci_table <- read.csv(paste0(root_dir,'outputs/multigraph_communities_small.csv'))

Next we subset the matrix by selecting the features under ‘miRNA category’. Each row in the Multi-Graph communities table table corresdponds to a feature annotation with non-zero assignments in the small round. The row index of each feature defines its column and row position in the Multi-Graph adjacency matrix.

## Extract microRNA matrix subset
btype <- 'miRNA'
idx <- which(multiloci_table$transcript_biotype == btype)
inM <- inputM[idx,idx]; ml <- multiloci_table[idx,]

Necessary steps to convert the adjacency matrix and the communities table to an igraph object are in the mg_build function.

## Generate microRNA graph object
attr <- 'transcript_name'
g = mg_build(inM, ml, attr)
g <- delete_vertices(g, V(g)[V(g)$weight < 1])

We can explore the patterns in the feature gene symbols established by the HGNC and link that to specific colors during visualization

## Extract microRNA mature extreme information from HUGO symol
V(g)$color_HUGO1 <- 'grey'
V(g)$color_HUGO1[grep('-3p',V(g)$feat)] <- 'orange'
V(g)$color_HUGO1[grep('-5p',V(g)$feat)] <- 'dodgerblue'

## Extract microRNA class from HUGO symbol
V(g)$color_HUGO2 <- 'grey'
idx <- grep('hsa-miR',as.character(V(g)$feat))
hugo <- as.factor(gsub('.*?([0-9]+).*', '\\1', V(g)$feat[idx]))
V(g)$color_HUGO2[idx] <- sample(grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = T)], 
                               length(unique(hugo)))[hugo]

The function mg_plotset generates a set of different plots of the Multi-Mapping graph and stores them as png files with a user-given file prefix (plotfile variable).

## Define plot colot
V(g)$customColor <- V(g)$color_HUGO1

## Standard plots
mg_plotset(g, plotfile= paste0(root_dir,'mgplots/mg_miRNA'))

Raw visualization of the graph

## Display plot
par(mar=c(0.01,0.01,0.01,0.01))
img <- magick::image_read(paste0(root_dir,'mgplots/mg_miRNA.png'))
plot(img)

Visualization of the graph colored by -3/-5 patterns in microRNA transcript symbol

## Display plot
par(mar=c(0.01,0.01,0.01,0.01))
img <- magick::image_read(paste0(root_dir,'mgplots/mg_miRNA_color.png'))
plot(img)

Visualization of the graph colored by -3/-5 patterns in microRNA transcript symbol and detected communities

## Display plot
par(mar=c(0.01,0.01,0.01,0.01))
img <- magick::image_read(paste0(root_dir,'mgplots/mg_miRNA_cl_color.png'))
plot(img)

mg_interactive(g, paste0(root_dir,'mgplots/mg_miRNA'))

We are done!